{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 基于脉冲的变分量子本征求解器\n",
    "\n",
    "\n",
    "*版权所有 (c) 2021 百度量子计算研究所，保留所有权利。*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 内容概要\n",
    "**注意：完整运行本教程的程序可能会花费超过 50 个 Quantum Hub 点数**\n",
    "\n",
    "本教程将介绍如何在脉冲层面实现变分量子本征求解器算法。本教程的大纲如下：\n",
    "\n",
    "- 变分量子本征求解器（VQE）\n",
    "- 基于脉冲的变分量子本征求解器（PBVQE）\n",
    "- 准备工作\n",
    "- 构造哈密顿量\n",
    "- 优化双量比特门的脉冲\n",
    "- 构造氢分子的哈密顿量\n",
    "- 构造基于脉冲的参数化电路及优化\n",
    "- 总结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 变分量子本征求解器（VQE）\n",
    "\n",
    "变分量子本征求解器（Variational Quantum Eigensolver, VQE）是在嘈杂中型量子（Noisy Intermediate-Scale Quantum, NISQ）计算机上运行的一种近似求解分子基态能量的算法。它的基本方法是估计给定哈密顿量的最小本征值并求其基态。对于近期的量子计算设备，门错误率较高、退相干时间较短以及连通性较差等问题限制了量子电路的深度。然而，VQE 算法只需要深度较低的量子电路即可实现，因而它被认为是利用 NISQ 设备解决实际问题的理想选择。\n",
    "\n",
    "VQE 的基本任务是制备参数化的量子态（trail state）$|\\psi(\\vec{\\theta})\\rangle$ 并估计出给定分子离散哈密顿量 $\\hat{H}_{\\rm mole}$ 的基态能量。其中，量子态 $|\\psi(\\vec{\\theta})\\rangle$ 是由参数化的量子电路（ansatz）生成。在这个过程中，我们采用经典的优化方法来寻找一组最优的参数 $\\vec{\\theta}^*$，以最小化期望值 $E = \\langle \\psi(\\vec{\\theta}) | \\hat{H}_{\\rm mole} | \\psi(\\vec{\\theta}) \\rangle$，即分子哈密顿量 $\\hat{H}_{\\rm mole}$ 的近似基态能量 $E_0^*$：\n",
    "\n",
    "$$\n",
    "E_0^* = {\\rm min}_{\\vec{\\theta}} \\langle \\psi(\\vec{\\theta}) | \\hat{H}_{\\rm mole} | \\psi(\\vec{\\theta}) \\rangle.\n",
    "$$\n",
    "\n",
    "在本教程中，我们将介绍在超导平台上使用 VQE 近似求解氢分子基态能量的基本方法。我们将考虑多种非理想因素，并从脉冲层面模拟 VQE 算法。首先，我们来介绍本教程中所使用的参数化量子电路模板，如下图所示：\n",
    "\n",
    "![VQE](figures/vqe-circuit-cn.png)\n",
    "\n",
    "它主要由若干参数化的单量子比特旋转门和 CNOT 门组成。由于 CNOT 门在超导平台上不能直接实现，因而在本教程中，我们将使用超导平台中实现效率更高的（hardware-efficient）双量子比特门，即 Cross-Resonance （CR）门来代替 CNOT 门作为纠缠门。同样，CR 门也可以配合若干单量子比特门产生最大纠缠态。理想 CR 门的矩阵为：\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\hat{U}_{\\rm CR}(\\alpha) = \\begin{bmatrix}\n",
    "\\cos{\\frac{\\alpha}{2}} & -i\\sin{\\frac{\\alpha}{2}} & 0 & 0 \\\\\n",
    "-i\\sin{\\frac{\\alpha}{2}} & \\cos{\\frac{\\alpha}{2}} & 0 & 0 \\\\ \n",
    "0 & 0 & \\cos{\\frac{\\alpha}{2}} & i\\sin{\\frac{\\alpha}{2}} \\\\\n",
    "0 & 0 & i\\sin{\\frac {\\alpha}{2}} & \\cos{\\frac{\\alpha}{2}} \n",
    "\\end{bmatrix}.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "在这里，我们设置 $\\alpha = -\\pi/2$。关于 CR 门的更多细节请[点击这里](https://quanlse.baidu.com/#/doc/tutorial-cr)。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 基于脉冲的变分量子本征求解器（PBVQE）\n",
    "\n",
    "在本教程中，我们将从脉冲层面研究 VQE 算法，我们称之为 **pulse-based VQE (PBVQE)**。与标准的 VQE 算法不同，PBVQE 不再优化逻辑量子电路中每个旋转门的参数，而是直接将脉冲参数作为优化参数来最小化损失函数（即基态能量）。下图显示了 PBVQE 和标准 VQE 算法之间的差异：\n",
    "\n",
    "![VQE](figures/vqe-scheme-cn.png)\n",
    "\n",
    "为了实现 PBVQE，我们需要将逻辑量子电路转换成**基于脉冲的参数化量子电路（pulse-based ansatz）**，即逻辑旋转门 $R_x(\\theta_m)$ 和 $R_y(\\theta_m)$ 分别被 $X$ 和 $Y$ 通道上振幅不同的控制脉冲所取代，我们称之为**基于脉冲的量子门（pulse-based gates）**：\n",
    "\n",
    "![VQE](figures/vqe-translate-cn.png)\n",
    "\n",
    "上图中，$U_{\\rm ENT}$ 是用于产生纠缠的酉算符（细节将会在后面的章节介绍）。这里，我们使用一种新的符号来表示**基于脉冲的量子门**的参数：\n",
    "\n",
    "$$\n",
    "\\vec{A} = [A_0, \\cdots, A_m, \\cdots, A_{M-1}],\n",
    "$$\n",
    "\n",
    "其中，$M$ 是**基于脉冲的量子门**的个数；$A_m$ 表示第 $m$ 个**基于脉冲的量子门**的高斯波形的振幅，因而脉冲包络的函数可以写为：\n",
    "\n",
    "$$\n",
    "\\Omega_m(t) = A_m e^{-(\\frac{t - \\tau_m}{\\sqrt{2} \\sigma_m}) ^2}.\n",
    "$$\n",
    "\n",
    "除脉冲强度以外的其它高斯脉冲参数，如宽度 $\\sigma_m$ 和中心位置 $\\tau_m$ 等在整个过程中都将被固定。这样一来，每个**基于脉冲的量子门**都只有一个需要优化的参数。引入**基于脉冲的参数化量子电路**后，在每次迭代中我们不再需要优化产生用于实现逻辑量子电路的驱动脉冲，这大大提高了 VQE 的效率和结果的准确性。\n",
    "\n",
    "在上面的章节中，我们简要介绍了传统的 VQE 和 PBVQE。在下面的部分中，我们将逐步演示使用量脉实现 PBVQE 的方法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 准备工作\n",
    "\n",
    "成功安装量脉后，您可以按照本教程运行下面的量脉程序。要运行此教程，您需要从量脉（Quanlse）和其它常用 Python 库导入以下包："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This module creates the Hamiltonian dictionary\n",
    "from Quanlse.QHamiltonian import QHamiltonian\n",
    "\n",
    "# These functions help us perform matrix calculations\n",
    "from Quanlse.Utils.Functions import tensor\n",
    "from Quanlse.Utils.Infidelity import unitaryInfidelity\n",
    "\n",
    "# These functions define useful operator matrices\n",
    "from Quanlse.QOperator import sigmaX, sigmaY, sigmaZ, sigmaI\n",
    "\n",
    "# This function generates wave data\n",
    "from Quanlse.QWaveform import QJob, QJobList, QWaveform, gaussian, square\n",
    "\n",
    "# This function uploads jobs to Quanlse Cloud Service and receives results\n",
    "from Quanlse.remoteSimulator import remoteSimulatorRunHamiltonian as runHamiltonian\n",
    "\n",
    "# This module defines matrices of the frequently used quantum gates\n",
    "from Quanlse.QOperation import FixedGate\n",
    "\n",
    "# This module saves the PBVQE results\n",
    "from Quanlse.Define import outputPath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the necessary packages\n",
    "import os\n",
    "from numpy import linalg, min, random, trace, dot, savez, load, identity, kron\n",
    "from math import pi\n",
    "from functools import reduce\n",
    "from scipy import optimize\n",
    "\n",
    "# Generate the path of npz file\n",
    "localFile = os.path.join(outputPath, f'pbvqe.npz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 构造哈密顿量\n",
    "\n",
    "首先，我们定义一些必要的常数，包括任意波形发生器（arbitrary wave generator, AWG）的采样周期、系统的量子比特的数量及能级。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sampling period (Nano second)\n",
    "dt = 2.0\n",
    "\n",
    "# Number of qubits\n",
    "qubits = 4\n",
    "\n",
    "# System energy level\n",
    "level = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后，我们定义超导量子比特的硬件参数。`freq` 列表中的项分别是 $\\omega_{\\rm q0}, \\omega_{\\rm q1}, \\omega_{\\rm q2}, \\omega_{\\rm q3}$，即每个量子比特的跃迁频率；`coupling` 列表中的项分别保存了量子比特 0-1、1-2、2-3、3-0 的耦合信息。利用旋转波近似（Rotating Wave Approximation, RWA），我们将系统定义在频率为 $\\omega_{\\rm RWA} = \\omega_{\\rm q0} = \\omega_{\\rm q2} = 4.914 \\times 2\\pi$ GHz 的旋转坐标系中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the hardware parameters of the qubits (GHz)\n",
    "freq = [4.914 * (2 * pi), 5.114 * (2 * pi), 4.914 * (2 * pi), 5.114 * (2 * pi)]\n",
    "\n",
    "# Define the coupling strength (GHz)\n",
    "coupling = [\n",
    "    [[0, 1], 0.016 * (2 * pi)],\n",
    "    [[1, 2], 0.016 * (2 * pi)],\n",
    "    [[2, 3], 0.016 * (2 * pi)],\n",
    "    [[3, 0], 0.016 * (2 * pi)]\n",
    "]\n",
    "\n",
    "# Frequency of rotating frame (GHz)\n",
    "rwa = 4.914 * (2 * pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们为所有的单量子比特门和双量子比特门设置固定的执行时间："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gate duration time (Nano second)\n",
    "tg2q = 200\n",
    "tg1q = 64"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随后，我们根据以下硬件结构使用量脉创建其哈密顿量，每个量子比特与其相邻量子比特耦合，耦合强度为一个恒定值：\n",
    "\n",
    "![VQE](figures/vqe-topo_structure-cn.png)\n",
    "\n",
    "上述系统的哈密顿量可以写为：\n",
    "$$\n",
    "\\hat{H}_{\\rm total} = \\sum_{q=0}^{3} \\delta_{q} \\hat{a}^{\\dagger}_{q}\\hat{a}_{q} + \\frac{1}{2}\\sum_{q=0}^{3}g_{q,(q+1) {\\rm\\ mod}\\ 4}(\\hat{a}_{q}\\hat{a}^{\\dagger}_{(q+1) {\\rm\\ mod}\\ 4}+\\hat{a}^{\\dagger}_{q}\\hat{a}_{(q+1) {\\rm\\ mod}\\ 4}) + \\sum_{q=0}^{3}\\Omega_{q}^x (t) \\hat{\\sigma}_{q}^{x} + \\sum_{q=0}^{3}\\Omega_{q}^y (t) \\hat{\\sigma}_{q}^{y} + \\sum_{q=0}^{3}\\Omega_{q}^z (t) \\hat{\\sigma}_{q}^{z} ,\n",
    "$$\n",
    "\n",
    "其中 $\\hat{a}_{q}$ 和 $\\hat{a}^{\\dagger}_{q}$ 分别是作用在第 $q$ 个量子比特的湮灭和产生算符。$\\hat{\\sigma}^x_{q}, \\hat{\\sigma}^y_{q}$ 和 $\\hat{\\sigma}^z_{q}$ 分别是作用在第 $q$ 个量子比特上的泡利算符。$\\delta_{q}=\\omega_{q} - \\omega_{\\rm RWA}$ 表示第 $q$ 个量子比特的失调强度；$g_{q,(q+1){\\rm\\ mod}\\ 4}$ 是第 $q$ 和第 $(q+1) {\\rm\\ mod}\\ 4$ 个量子比特之间的耦合强度； $\\Omega_q^{x,y,z}(t)$ 是作用在第 $q$ 个量子比特上的磁通调控或微波调控的包络函数。我们可以使用量脉方便地定义上述系统的哈密顿量："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the Hamiltonian\n",
    "vqeHam = QHamiltonian(qubits, level, dt)\n",
    "\n",
    "# Add the coupling terms\n",
    "for item in coupling:\n",
    "    q0, q1 = item[0][0], item[0][1]\n",
    "    vqeHam.addCoupling([q0, q1], g=item[1] / 2)\n",
    "\n",
    "for qubit in range(qubits):\n",
    "    # Add the detuning terms\n",
    "    detuning = freq[qubit] - rwa\n",
    "    vqeHam.addDrift(sigmaZ, qubit, coef=detuning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "关于使用量脉构建哈密顿量的更多方法，可以查看教程[单量子比特门](https://quanlse.baidu.com/#/doc/tutorial-single-qubit)。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 优化双量子比特门\n",
    "\n",
    "在本教程中，我们使用 CR 门作为纠缠门（关于 CR 门的更多信息，可以查看教程：[Cross-Resonance 门](https://quanlse.baidu.com/#/doc/tutorial-cr)）。由于在本教程中，相邻量子比特之间的耦合方式为直接耦合，因此在一个量子比特上施加 $X$ 脉冲会同时影响到另外两个相邻的量子比特。因此，我们在设计脉冲时需要考虑这个因素，以抑制串扰造成的影响。\n",
    "\n",
    "![VQE](figures/vqe-crosstalk-cn.png)\n",
    "\n",
    "在这里，我们使用 `vqeHam.subSystem()` 从系统哈密顿量 `vqeHam` 中提取两个由三个量子比特组成的子系统用于优化 CR 门，其中一个是由量子比特 0-1-2 组成的子系统，另一个是由量子比特 1-2-3 组成的子系统。在这些子系统上，我们分别设置 $\\hat{U}_{\\rm goal}=I\\otimes\\hat{U}_{\\rm CR}$ 作为目标酉矩阵来优化相应脉冲，即在子系统的第二和第三个量子系统上生成一个 CR 门。\n",
    "\n",
    "我们定义函数 `makeCrPulse()` 用于生成 CR 门所需的脉冲序列。我们在当前子系统的第二个量子比特上施加高斯微波驱动脉冲，同时固定其宽度和中心位置，并将其振幅作为优化的第一个参数。第二个需要优化的参数是施加在第一个量子比特上的磁通控制的振幅。请注意，`tag=\"det\"` 的驱动还同时用于将旋转坐标参考系转换为特定频率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeCrPulse(ham, subSys3q, driveFreq, amp, shift, t):\n",
    "    \"\"\" Assemble the pulses for CR gates \"\"\"\n",
    "    subHam = ham if subSys3q is None else ham.subSystem(subSys3q)\n",
    "    subHam.clearWaves()\n",
    "    subHam.appendWave(sigmaX, 1, gaussian(t, amp, tg2q / 2, tg2q / 8), tag=\"XY\")\n",
    "    # frame transformation\n",
    "    subHam.appendWave(sigmaZ, 0, square(t, rwa - driveFreq + shift), tag=\"Z\")\n",
    "    subHam.appendWave(sigmaZ, 1, square(t, rwa - driveFreq), tag=\"det\")\n",
    "    subHam.appendWave(sigmaZ, 2, square(t, rwa - driveFreq), tag=\"det\")\n",
    "    return subHam.job if subSys3q is None else subHam.outputInverseJob(qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随后，我们定义一个函数 `optimize_cr()` 来进行优化过程，并保存最佳参数以供进一步使用。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimizeCr(subSys3q, driveFreq):\n",
    "    \"\"\" Realize a CR gate on the second & third qubits \"\"\"\n",
    "    crHam = vqeHam.subSystem(subSys3q)\n",
    "    uGoal = tensor([identity(2), FixedGate.CR.getMatrix()])\n",
    "\n",
    "    def crLoss(_x):\n",
    "        # Clear and add waves\n",
    "        crHam.clearWaves()\n",
    "        # Generate and add waves for CR gate implementation\n",
    "        _crJob = makeCrPulse(crHam, None, driveFreq, _x[0], _x[1], tg2q)\n",
    "        # Simulate the system's evolution and obtain the infidelity\n",
    "        unitary = crHam.simulate(job=_crJob)[0][\"unitary\"]\n",
    "        infidelity = unitaryInfidelity(uGoal, unitary, 3)\n",
    "        return infidelity\n",
    "\n",
    "    opt = optimize.dual_annealing(crLoss, [(-2, 2), (-0.2, 0.2)], maxiter=60)\n",
    "    print(\"Min infidelity:\", opt[\"fun\"])\n",
    "    return opt[\"x\"][0], opt[\"x\"][1]\n",
    "\n",
    "lhlQ1X, lhlQ0Z = optimizeCr([0, 1, 2], 4.914 * 2 * pi)\n",
    "hlhQ1X, hlhQ0Z = optimizeCr([1, 2, 3], 5.114 * 2 * pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 构造氢分子的哈密顿量\n",
    "\n",
    "在这一节中，我们将介绍如何在脉冲层面上估计氢分子的基态能量。我们将省略费米子—量子比特（fermion-to-qubit）映射的具体细节（请访问[量桨](https://github.com/PaddlePaddle/Quantum/blob/master/tutorial/quantum_simulation/VQE_CN.ipynb)获得更多相关信息）。首先，我们定义一个函数 `pauli_str_to_matrix()`，将**泡利字符串**转换为氢分子的离散哈密顿量 $\\hat{H}_{\\rm mole}$："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauliStrToMatrix(pauli_str, n):\n",
    "    \"\"\"\n",
    "    Convert the Pauli string in Hamiltonian\n",
    "    \"\"\"\n",
    "    def NKron(AMatrix, BMatrix, *args):\n",
    "        return reduce(\n",
    "            lambda result, index: kron(result, index),\n",
    "            args,\n",
    "            kron(AMatrix, BMatrix), )\n",
    "    pauli_dict = {\n",
    "        'i': sigmaI().matrix,\n",
    "        'x': sigmaX().matrix,\n",
    "        'y': sigmaY().matrix,\n",
    "        'z': sigmaZ().matrix\n",
    "    }\n",
    "    # Parse pauli_str; 'x0,z1,y4' to 'xziiy'\n",
    "    new_pauli_str = []\n",
    "    for coeff, op_str in pauli_str:\n",
    "        init = list('i' * n)\n",
    "        op_list = op_str.split(',')\n",
    "        for op in op_list:\n",
    "            pos = int(op[1:])\n",
    "            assert pos < n, 'n is too small'\n",
    "            init[pos] = op[0]\n",
    "        new_pauli_str.append([coeff, ''.join(init)])\n",
    "\n",
    "    # Convert new_pauli_str to matrix; 'xziiy' to NKron(x, z, i, i, y)\n",
    "    matrices = []\n",
    "    for coeff, op_str in new_pauli_str:\n",
    "        sub_matrices = []\n",
    "        for op in op_str:\n",
    "            sub_matrices.append(pauli_dict[op])\n",
    "        if len(op_str) == 1:\n",
    "            matrices.append(coeff * sub_matrices[0])\n",
    "        else:\n",
    "            matrices.append(coeff * NKron(sub_matrices[0], sub_matrices[1], *sub_matrices[2:]))\n",
    "\n",
    "    return sum(matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这里，我们使用原子间隔为 $d=74$ pm 的氢分子空间构型数据，这些数据来自[量桨](https://github.com/PaddlePaddle/Quantum/blob/master/tutorial/quantum_simulation/VQE_CN.ipynb)。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targetHam = [\n",
    "    [-0.042078976477822, 'i0'],\n",
    "    [ 0.177712874651399, 'z0'],\n",
    "    [ 0.177712874651399, 'z1'],\n",
    "    [-0.242742805131446, 'z2'],\n",
    "    [-0.242742805131462, 'z3'],\n",
    "    [ 0.170597383288005, 'z0,z1'],\n",
    "    [ 0.044750144015351, 'y0,x1,x2,y3'],\n",
    "    [-0.044750144015351, 'y0,y1,x2,x3'],\n",
    "    [-0.044750144015351, 'x0,x1,y2,y3'],\n",
    "    [ 0.044750144015351, 'x0,y1,y2,x3'],\n",
    "    [ 0.122933050561837, 'z0,z2'],\n",
    "    [ 0.167683194577189, 'z0,z3'],\n",
    "    [ 0.167683194577189, 'z1,z2'],\n",
    "    [ 0.122933050561837, 'z1,z3'],\n",
    "    [ 0.176276408043195, 'z2,z3']\n",
    "]\n",
    "hMatrix = pauliStrToMatrix(targetHam, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述分子哈密顿量基态能量的理论值可以通过如下方法计算："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the theoretical eigenvalue\n",
    "eigVal, eigState = linalg.eig(hMatrix)\n",
    "minEigH = min(eigVal.real)\n",
    "print(f\"Ground state energy: {minEigH} Ha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 基于脉冲的量子电路\n",
    "\n",
    "首先，我们参考标准 VQE 中最常用的参数化量子电路模板，设计了一个基于脉冲的量子电路。下图显示了该量子电路中的一层，其中，每个量子比特都包含 3 个单量子比特门，而每个单量子比特门都有一个参数作为高斯脉冲包络的最大振幅，脉冲宽度和中心位置是固定的。\n",
    "\n",
    "![VQE](figures/vqe-scheduling-cn.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于脉冲电路较为复杂，因而我们定义了一个函数 `makeWaveSchedule()` 专门用于生成并排列上述电路所对应的脉冲序列。其中，参数 `x` 是优化参数列表（即脉冲参数 $\\vec{A}$）；`vqeJob` 是由 `addWave()` 生成的波形数据列表，用于保存用户定义波形的详细信息。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeWaveSchedule(x):\n",
    "    \"\"\" Generate waves for pulse-based circuit \"\"\"\n",
    "    # Generate pulses for CR gate\n",
    "    crJob = vqeHam.createJob()\n",
    "    crJob += makeCrPulse(vqeHam, [3, 0, 1], 5.114 * 2 * pi, hlhQ1X, hlhQ0Z, tg2q)\n",
    "    crJob += makeCrPulse(vqeHam, [0, 1, 2], 4.914 * 2 * pi, lhlQ1X, lhlQ0Z, tg2q)\n",
    "    crJob += makeCrPulse(vqeHam, [1, 2, 3], 5.114 * 2 * pi, hlhQ1X, hlhQ0Z, tg2q)\n",
    "    crJob += makeCrPulse(vqeHam, [2, 3, 0], 4.914 * 2 * pi, lhlQ1X, lhlQ0Z, tg2q)\n",
    "    # Assemble the pulses\n",
    "    depth = int(len(x) / 12)\n",
    "    vqeJob = vqeHam.createJob()\n",
    "    for d in range(depth):\n",
    "        gate1QJob = vqeHam.createJob()\n",
    "        # Add pulses for single-qubit gates\n",
    "        for q in range(4):\n",
    "            # X/Y/X controls\n",
    "            gate1QJob.addWave(sigmaX, q, gaussian(tg1q, x[12 * d + q], tg1q / 2, tg1q / 8), t0=0)\n",
    "            gate1QJob.addWave(sigmaY, q, gaussian(tg1q, x[12 * d + 4 + q], tg1q / 2, tg1q / 8), t0=tg1q)\n",
    "            gate1QJob.addWave(sigmaX, q, gaussian(tg1q, x[12 * d + 8 + q], tg1q / 2, tg1q / 8), t0=tg1q * 2)\n",
    "            # Set detuning\n",
    "            gate1QJob.addWave(sigmaZ, q, square(tg1q * 3, rwa - freq[q]), t0=0, tag=\"det\")\n",
    "        vqeJob += gate1QJob\n",
    "        vqeJob += crJob\n",
    "    return vqeJob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本教程中，我们使用 `Scipy` 提供的基于梯度的优化方法（L-BFGS-B）来最小化目标函数。在每次迭代中，L-BFGS-B 需要用户提供每个参数的梯度信息，在这里我们使用两点有限差分法来近似计算梯度：\n",
    "$$\n",
    "\\frac{\\partial{\\rm Loss}(\\vec{A})}{\\partial a_m} = \\frac{{\\rm Loss}(a_0, \\cdots, a_m + \\epsilon, \\cdots, a_{M-1}) - {\\rm Loss}(a_0, \\cdots, a_m - \\epsilon, \\cdots, a_{M-1})}{2\\epsilon} ,\n",
    "$$\n",
    "\n",
    "其中，$\\vec{A} = [A_0, \\cdots, A_{M-1}]$ 是脉冲参数列表，$\\epsilon$ 是一个很小的正数，而损失函数 ${\\rm Loss}(\\vec{A})$ 定义为：\n",
    "\n",
    "$$\n",
    "{\\rm Loss}(\\vec{A}) =  \\langle \\psi(\\vec{A}) | \\hat{H}_{\\rm mole} | \\psi(\\vec{A}) \\rangle.\n",
    "$$\n",
    "\n",
    "其中，量子态 $\\psi(\\vec{A})$ 是基于脉冲的量子电路所产生的。有限差分法需要大量的样本，例如，当脉冲参数的参数为 $M$ 时，我们需要 $2M$ 次采样来估计近似梯度。因此，我们使用量脉云服务来加速这个过程。\n",
    "\n",
    "为了使用量脉云服务，我们需要导入 `Define` 并传入 token，用户可以在 [Quantum-hub](http://quantum-hub.baidu.com) 申请获得 token。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the loss function\n",
    "import copy\n",
    "from Quanlse import Define\n",
    "Define.hubToken = \"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们定义 VQE 的损失函数。在这个函数中，我们模拟了脉冲参数为 $\\vec{x}$ 时基于脉冲的电路的演化，并用上面提到的有限差分法计算这一点的梯度。在每次迭代中，我们将当前脉冲参数列表 $\\vec{x}$ 输入到损失函数中，并将所有采样所需的脉冲数据生成并打包到 `waveList` 中。最终，`waveList` 包含用于求解梯度的 $2M$ 次采样和用于获取损失值的 1 个采样的脉冲数据。\n",
    "\n",
    "我们在上面的步骤中将所有的任务集成到一个列表中，即 `waveList`，并通过函数 `runHamiltonian()` 将任务列表提交给量脉云服务。正常情况下，大约 15 到 20 秒后，我们将收到返回结果，结果将作为 JSON 文件保存到 `Output` 文件夹中。同时，变量 `result` 会被赋予一个列表，其中包含与 `waveList` 对应的所有模拟结果。\n",
    "\n",
    "**注意**：`waveList` 的每一项都包含由 `makeWaveSchedule()` 函数生成的基于脉冲的 VQE 的所有波。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss(x):\n",
    "    global lossHistory\n",
    "    # Add wave for current point\n",
    "    waveList = vqeHam.createJobList()\n",
    "    waveList.addJob(makeWaveSchedule(x))\n",
    "\n",
    "    # Add wave for calculating gradient\n",
    "    for xId in range(len(x)):\n",
    "        xList = copy.deepcopy(x)\n",
    "        xList[xId] -= 1e-8\n",
    "        waveList.addJob(makeWaveSchedule(xList))\n",
    "        xList[xId] += 2 * 1e-8\n",
    "        waveList.addJob(makeWaveSchedule(xList))\n",
    "\n",
    "    # Simulate the evolution\n",
    "    result = runHamiltonian(vqeHam, jobList=waveList)\n",
    "\n",
    "    # Calculate the loss function\n",
    "    lossList = []\n",
    "    for item in result:\n",
    "        state = item[\"unitary\"]\n",
    "        lossVal = (state.conj().T @ hMatrix @ state).real[0][0]\n",
    "        lossList.append(lossVal)\n",
    "\n",
    "    # Calculate the gradients\n",
    "    gradient = []\n",
    "    for index in range(len(x)):\n",
    "        gradient.append((lossList[2 + 2 * index] - lossList[1 + 2 * index]) / 1e-8 / 2)\n",
    "\n",
    "    print(\"Loss function:\", lossList[0])\n",
    "    lossHistory.append(lossList[0])\n",
    "    return lossList[0], gradient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后我们使用由 `Scipy` 提供的 `fmin_l_bfgs_b()` 函数最小化前面定义的损耗函数。\n",
    "\n",
    "**注意**：此优化可能需要超过 15 分钟。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth = 3\n",
    "lossHistory = []\n",
    "initParas = [random.rand() for _ in range(depth * 12)]\n",
    "bounds = [(-1.5, 1.5) for _ in range(depth * 12)]\n",
    "x, f, d = optimize.fmin_l_bfgs_b(loss, initParas, fprime=None, bounds=bounds, maxiter=200)\n",
    "\n",
    "# Save the loss history to a file for further usage\n",
    "savez(localFile, lossHistory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"The estimated ground state energy is: {f} Ha\")\n",
    "print(\"Total iteration:\", d[\"nit\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可见，最终收敛的精度很高，迭代次数为 72 次。随后，我们绘制完整的迭代过过程："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the loss_history list from the npz file.\n",
    "lossHistory = load(localFile)['arr_0']\n",
    "\n",
    "# Plot the figures\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(range(len(lossHistory)), lossHistory, label=\"Energy\")\n",
    "plt.axhline(minEigH, c=\"gray\", ls=\"--\", lw=1.0)\n",
    "plt.xlabel(\"Iteration\")\n",
    "plt.ylabel(\"Energy (Ha)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，我们可以使用 `plot()` 方法绘制脉冲序列："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the waveforms.\n",
    "makeWaveSchedule(x).plot(color=['red', 'green', 'blue'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 总结\n",
    "用户可以通过点击这个链接 [tutorial-pbvqe.ipynb](https://github.com/baidu/Quanlse/blob/main/Tutorial/CN/tutorial-pbvqe-cn.ipynb) 跳转到此 Jupyter Notebook 文档相应的 GitHub 页面并获取相关代码以运行该程序。我们鼓励用户使用量脉开发更多脉冲层的 NISQ 算法。\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "\\[1\\] [Peruzzo, Alberto, et al. \"A variational eigenvalue solver on a photonic quantum processor.\" *Nature communications* 5 (2014): 4213.](https://doi.org/10.1038/ncomms5213)\n",
    "\n",
    "\\[2\\] [Moll, Nikolaj, et al. \"Quantum optimization using variational algorithms on near-term quantum devices.\" *Quantum Science and Technology* 3.3 (2018): 030503.](https://doi.org/10.1088/2058-9565/aab822)\n",
    "\n",
    "\\[3\\] [Kandala, Abhinav, et al. \"Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.\" *Nature* 549.7671 (2017): 242-246.](https://doi.org/10.1038/nature23879)\n",
    "\n",
    "\\[4\\] [Rigetti, Chad, and Michel Devoret. \"Fully microwave-tunable universal gates in superconducting qubits with linear couplings and fixed transition frequencies.\" *Physical Review B* 81.13 (2010): 134507.](https://doi.org/10.1103/PhysRevB.81.134507)\n",
    "\n",
    "\\[5\\] [Meitei, Oinam Romesh, et al. \"Gate-free state preparation for fast variational quantum eigensolver simulations: ctrl-VQE.\" *arXiv preprint arXiv:2008.04302* (2020).](https://arxiv.org/abs/2008.04302)\n",
    "\n",
    "\\[6\\] [Wilhelm, Frank K., et al. \"An introduction into optimal control for quantum technologies.\" *arXiv preprint arXiv:2003.10132* (2020).](https://arxiv.org/abs/2003.10132)\n",
    "\n",
    "\\[7\\] [Krantz, Philip, et al. \"A quantum engineer's guide to superconducting qubits.\" *Applied Physics Reviews* 6.2 (2019): 021318.](https://aip.scitation.org/doi/abs/10.1063/1.5089550)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
